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SECTION I 
INTRODUCTION 


Astronomy has traditionally made great use of imagery, but this 
data has almost exclusively been evaluated via the photographic process. 
With the advent of the unmanned planetary missions and digital cameras, 
techniques for processing digital Images have been developed and are 
beginning to be applied in astronomy. The advantages of digital images 
over photographic ones are very great and, with the development of new 
algorithms and cameras, become greater all the time. It is thus of im- 
portance to search for applications and disseminate digital techniques 
to the astronomical community at large. When Space Telescope becomes 
operational, the quantity and unprecedented quality of the calibrated 
digital images returned will call for every available processing means 
at our disposal. This report is intended to provide a few of those 
means. 

This work was performed during 1979 at the Image Processing Labora- 
tory (IPL) of JPL under FY79 RTOP #389— M 1—0 1 . It represents a continuing 
effort (References 1 and 2) to identify and develop new image processing 
applications to astronomy. Additional goals of this research are the 
further education of the astronomical community in digital image processing 
techniques and the expansion of methods for analysis of high-quality space 
telescope imagery when it becomes available. Hopefully, others will be 
inspired to find applications to their own disciplines. 

This report is divided into eight sections. Sections II through VII 
present techniques applied to a major topic. The data used and the tech- 
niques involved are described in detail, along with photographs and plots 
of the results. Section VIII gives conclusions that can be made as a 
result of this work. This report illustrates image processirg and does 
not interpret the physical significance of results that develop from the 
data. The role of data Interpretation resides with the collaborating 
astronomer. Image processing is a vehicle to strengthen and simplify 
his task. 


1-1 




SECTION II 


DECONVOLUTION - MAXIMUM ENTROPY VS. WIENER 


The algorithm developed by Dr. Gull end Dr. Deniell (References 
3 end 4) for eorreoting incomplete radio naps was provided in a fora 
suitable for deconvolution of optioal inages (Reference 5). This tech- 
nique has proven very successful in restoring Fourier transforms, and 
the philosophy is applicable to solving, in an iterative fashion, the 
general inverse problem. 

The maximum entropy technique has several desirable properties 
which are not observed by linear filters: 

(1) It produces a unique map. 

(2) The restored intensity values are non-negative. 

(3) Ringing (Gibbs phenomena) is virtually undetectable. 

To provide a control, an 8-bit synthetic image was created (upper 
left corner of Figure 1). This image was convolved with a 5-by-5 element 
kernel of equal weight which replaced each pixel by the means of its 
neighborhood (upper right of Figure 1). The objeot was to restore 
the degraded image to its original form. In the central row of Figure 1, 
the left image represents the maximum entropy restoration and the right 
image represents the Wiener linear restoration (References 6 and 7). 

The lower left image of Figure 1 contains the chi squared residuals 
map of the maximum entropy restoration, i.e., the ratio of misfit to 
the data to the error due to noise in the data. The lower right image 
of Figure 1 contains the log MTF of the kernel used to degrade the 
original, i.e., the log modulus of the Fourier transform of the 5-by-5 
element point spread function. T' t a particular kernel is not trivial 
to correct for because it introduc s both 180° phase shifts and crosses 
the zero amplitude plane (dark bars) in many places. 

In the case illustrated by Figure 1 (which is noiseless if one 
can disregard the RMS integer truncation noise of 1 /2<fl ), the Wiener 
model, even with a modest signal -to -noise ratio of 15, is more effective 
for visual recognition. 

Figures 2 and 3 demonstrate the effect of adding random noise 
to the degraded image before restoration in RMS quantities of 30 and 
60, respectively. In Figure 3, the true signal-to-noise ratio in the 
image is approximately unity and the linear model has been forced to 
further blur the data in order to improve it. The two algorithms operate 
with equal effectiveness in Figure 2. In Figure 3, the noise level is 
almost too high for either algorithm to effect an improvement. 

In Figure 4, a larger filter of 9-by-9 elements has been employed, 
which has destroyed the control beyond the human pattern recognition 
capability to recognize any characters. No noise was added and the 
linear model was able to restore most of the image into the recognizable 
range. 
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It appears fro® th« above exercise that the naxiau® entropy resto- 
ration is able to coopensate for snail amounts of blurring in the presei.o 
of noise and to produoe relatively artifact- free inages when the noise 
level is moderate. The Wiener model is superior fro® a subjective 
standpoint only in the low-noise regime. 

To compare both algorithms on astronomical images, four plates 
taken through a Polaroid at P.A. 245° of M87 were averaged. Reference 
stars were selected from the field and placed into a row at the bottom 
of the image. These stars provided the required point spread function. 
Figure 5 illustrates a mosaic of several Images of M87. On the top 
row the left image is the original, the middle is the maximum entropy 
restored version, and the right image is the chi squared residuals 
map. On the bottom row from left to right are linear Wiener model 
restorations with slgnal-to-noiae values of 4, 8, and 14. The maximum 
entropy restoration is virtually artifact free but is only slightly 
enhanced. The Wiener restorations, although full of artifacts, provide 
better restoration where the signal- to-noiae ratio is high. The back- 
ground artifact is due to the linear model's assumption that the 
signal-to-noise ratio is independent of position. It is ironic that 
the chi squared residuals map is as interesting as the restoration. 

In Figure 6, a similar picture to Figure 5 is presented, which provides 
a better maximum entropy restoration. This was accomplished by smoothing 
the original image with a small convolutional filter, thus improving 
the overall signal- to-noise ratio by suppressing the highest frequencies. 

It is important to consider that, because the maximum entropy 
technique produces little artifact, it may be more valuable to employ 
in cases where the true nature of the object cannot be estimated beforehand. 

One of the powerful capabilities of the maximum entropy restoration 
is its ability to assign noise weighting to the image on a pixel- to- pixel 
basis. As an example of this, the control image used in Figure 1 was 
degraded with a 7-by-7 kernel and then decimated with a random grid (shown 
at the bottom left of Figure 7). The lower right of Figure 7 shows 
the restoration. Although this restoration is not better than the blurred 
version, the application of any linear convolutional model would have 
resulted in a totally unrecognizable picture. 

A comparison of both algorithms was also performed upon synthetic 
spectra in order to evalute their performance in one dimension. Each 
of the three figures, Figures 6, 9, and 10, contains the original spectrum 
at the top, the degraded and noise added spectrum at the bottom and, 
in the second and third rows from the top, the maximum entropy and Wiener 
model restorations respectively. The object was, as usual, to reconstruct 
the top row from the bottom one. Figures 8, 9, and 10 differ insofar 
as the KM3 noise added to the degraded spectrum was zero for Figure 8; 
two for Figure 9; and six for Figure 10; with corresponding Wiener 
signal-to-noise values of 20, 10, and 4. From a comparison of Figures 8, 

9, and 10, it can be seen that once again as in the case of the two- 
dimensional restorations the linear model is superior when the noise 
level is very low. At intermediate noise levels, the two algorithms 
are very similar, and, at high noise levels, they become equally unreli- 
able, producing artifacts. 
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Figure 2. Maximum entropy vs. Wiener restoration on a blurred target 

with an RMS = 30 noise level .idded. The Wiener signal-to-noise 
ratio was 5. 
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Figure 6. Maximum entropy vs. Wiener restoration on a real scene of 

M87 : (top left) smoothed original; (middle) maximum entropy 

restoration; (right) chi squared residuals. 
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Maximum entropy restoration of blurred decimated images: 

(top left) original; (top right) blurred image; (lower left) 
decimated image; (lower right) restoration. 





Figure 8. Maximum entropy vs. Wiener restorations for one-dimensional 

spectra: (top) original spectrum; (bottom) noiseless degraded 

spectrum; (second row) maximum entropy restoration; (third 
row) Wiener restoration with a signal-to-noise ratio of 20. 
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Figure 10. Maximum entropy vs. Wiener restoration for an RMS noise 

of 6 added to the degraded image. The Wiener signal-to-noise 
ratio used was 4. 
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SECTION III 


POLARIZATION IN GALAXIES 


Plates of M82 obtained by Dr. Allan Sandage of Hale Observatories 
in two polarization orientations were scaled such that the nean sky 
level was centered at_zero and the variances of both images were identical. 
Even though only two 3-vector orientations were known, it was possible 
to compute the polarization beoause the direction of the polarization 
vector was known to be approximately tangential about the center of 
the M82 nucleus (Reference 8}. 

Figure 11 illustrates the difference between two images in the 
sense 3-vector along major axis minus 3-vector along minor axis, feich 
filamentary structure is visible, all due to polarization. For the 
case of four 3-vector positions 45° apart, the polarization can be 
computed from Equations 1 through 4. 


G 

I 


I 0° “ *90° 
x 0° + I 90° 


(D 


U x 45° “ I 135° 

mm JJ 

I Ij, 5 o + I 135 o 


tan 20 s 


U 

Q 


P 



(3) 


(4) 


Relations (Eq. 1-;) can be combined to give 


J tan 26 


U 

I 


(5) 


and from this and relation (Eq. 4) 


P 


Q 

i 


Y 1 + tan 2 20 


( 6 ) 


3-1 


the form when the angle is known. The tan 2 20 acts as a modulating term 
which corrects for the polarization amplitude at other that 0 and 90 
degrees. At multiples of 90° beginning at 45°, the function becomes 
unstable and results in noise amplification, leaving a cross in the 
polarization picture. Figures 12 and 13 illustrate the polarization 
map of M82 both before and after local smoothing of the original Images. 
In both cases 25 percent polarization is equivalent to white, and zero 
is equivalent to black. Because of the low signal-to-noise ratio near 
the plate background, low-intensity values become corrupted by noise 
and result in spuriously high polarization values giving rise to the 
white ring around M82. A threshold was applied to the algorithm which 
set the polarization to zero if the Intensity in either frame fell below 
a certain level. This gave rise to the dark region outside the bright 
noise-induced halo. Due to assumptions made in scaling the images, 
the polarization on the bar of the galaxy is constrained to be zero. 

The filaments, which extend normal to the bar, are highly polarized 
at about 20 percent. 

From inspection of the previous formula, it is clear that noise 
will severely influence the measure of polarization if either: 

(1) The differences between intensities are on the same order 
as the RMS noise. 

(2) The intensities are so low that their sums are on the same 
order as the RMS noise. 



Figure 11 . M82 difference image in the sense E-vector along the major 

axis minus E-vector normal to the major axis 
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hi^ure it’. 'Ao2 polarization image scaled with 2b percent polarization 

to white. The cross is artifact due to the assumed Circumferen- 
tial direction of the polarization vector. The unite halo around 
Hbi; is artifact due to low signal-to-noise ratio close 
to tne swy . Polarized filaments are visible in two quadrants 
above and below tne cross center. 
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Figure 13* M62 polarization image in the same format as Figure 12 , 

but computed from .smoothed images 
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SECTION IV 


TIME CHANGES IN M87 


Plat** of M87 were taken In two epochs, on* in 1956 by Dr. Baad* 
and tha othar in 1978 by Dra. Arp and Sulantio of Halo Obsarvatoriaa. 

Eaoh aat of plataa waa racordad at 3* and 6-minute axpoauraa and aaoh 
at four polarisation oriantationa. A mosaic of all sixtaan images 
is praaantad in Pigura 1*. Eaoh image waa aubjaotad to an approximate 
HD oorraotion by foroing tha atallar profila* to ba Gaussian, thara 
baing no othar calibration available. Eaoh of tha image* was than 
ragistarad using foreground stars, and several roforanca stars ware 
collected into a row below tha corresponding image of M87. 

In order to observe any relative intensity variations with time 
in tha knots along tha main jet, tha images war* normalized such that 
tha oentral galaxy and its halo war* of the same amplitude. A similar 
analysis but normalizing tha rafaranca stars to the asm* amplitude 
was performed and is in press (Reference 9). After normalizing tha 
halo, diffaranoa pictures in tha sense of 1978 compared to 1956 ware 
genera lad and are presented as Figures 15 and 16; Figure 15 is tha 
difference between tha 6-minute exposures, and Figure <6 is tha difference 
between tha 3-alnute exposures. It is clear from tha diffaranoa pictures 
that at least the inner main knot "A" has increased in brightness relative 
to the halo from 1956 to 1978, especially at polarization P.A. 335°. 




Figure 14. M87 mosaic displaying four polarization orientations, each 

in 6- and 3-minute exposures and each at two epochs, 1956 
and 1978 
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Figure lb. M87 difference images in the sense 1^78 minus 'dSb, each 

in four polarization orient at ions . The images were scaled 
so that the halo of M 87 matched in intensity. Heferenoe 
stars appear below. These are b-minute exposures. 
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Figure 16. M87 difference images in the same format as Figure 15. 

Thesp are 3-mirute exposures. 
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SECTION V 


PLANETARY NEBULAE - COMPARING EMISSION LINE IMAGES 


In this section, the steps leading toward calculation of the 
volume emission coefficient, temperature, and electron density from 
emission line images of planetary nebulae are described. 

The work in this section was inspired and performed by Peter 
Kupferman of JPL, using the facilities of the Image Processing Laboratory. 
Images were obtained at the Palomar 60-inch telescope with a silicon 
vidicon through interference filters, listed in Table 1, of three 
planetarie3 NGC7662, NGC6720, and NGC2392. The images comprised 
256-by-256 pixels of 12 bits each at 0.5 seconds of arc per pixel. 

Each image was radiometrically decalibrated by subtracting averaged 
dark current frames and dividing by averaged flat fields obtained during 
the observing run, implying a linear light transfer curve model. All 
of the images for each object were then geometrically registered using 
nearby reference stars unrelated to the nebula. All of the images 
obtained at each wavelength (there were from 2 to 9 redundant pictures) 
were averaged to provide a single picture at each wavelength. A low 
spatial frequency analytical model was generated to match the residual 
background sky level after decalibration. This model was subtracted 
from each image in order to reduce small inhomogeneities in the sky. 

Figure 17 illustrates a mosaic of each of the images for NGC6720 at 
this stage in the processing. Each image is arbitrarily scaled to 
produce a visible pictur . Figure 18 displays the log amplitude of the 
Fourier transform of each of the images in Figure 17. It is evident 


Table 1. Wavelengths of Each Planetary Nebula Image 


Target Line I.D. 

A 

0 

Bandpass A 

ton] 

3727 

20 

[OIII] 

4363 

13 

Hell 

4686 

28 

H(3 

4861 

15 

Hel 

5876 

15 

[OIII] 

5007 

17 

[01] 

6300 

6 

Ha 

6563 

11 

Ha+[NII] 

6563 

175 

[SII] 

6716 

10 

[SII] 

6731 

11 
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from these figures that a periodic pattern due to the deceleration 
grid in the vidicon (hash lines) and clusters due to quasi-coherent 
noise are still present to a varying degree in each image. From 
inspection of Figure 18, it was determined that the system MTF went 
to zero at about 0.25 cycles per sample due to atmospheric seeing; 
thus, it was possible to reduce the coherent, periodic, and random noise 
greatly by convolving with a band pass filter that suppressed information 
above this frequency. This technique obviated the need for a custom 
modification of each Fourier transform. Figure 19 displays each of 
the rescaled Fourier transforms after convolution, and Figure 20 
illustrates the corresponding smoothed images. 

In order to compute the volume emission coefficient, each image 
was first converted into polar coordinates (Figure 21) in which each 
row represents a radial section outward from the central star, the 
rows extending in azimuth from 0° at the top to 360° at the bottom, 
if one can assume spherical symmetry, then each row in the polar 
coordinate space would provide a separate model for the nebula, allowing 
different models for different structures at those azimuths. From 
the observed projected radial intensity profiles I(x), one can obtain 
the true radial intensity profile j(r), known as the volume emissivity 
from the Abel transform: 



From the true radial emission j(r), one can regenerate the nebula by 
using the inverse polar coordinate transformation. In Figure 22, the 
positive volume emission coefficient is displayed in red with negative 
regions coded in green. Negative emission indicates where the planetary 
.. ^del is not circum ferentially symmetric and that a more sophisticated 
model is desired. 

Absolute calibrations were obtained from spectrophotometry performed 
at specific points in the nebula at the McGraw Hill 52-inch telescope 
at Kitt Peak. The intensities at these points could be used to scale 
tne digital images linearly so that each pixel was known in units of 
ergs/cm^-steradiar-sec , providing units of j(r) in ergs/cm^-sterad ian-sec . 
rtelative photometry is good from 2 to 4 percent, and absolute photometry is 
good from 5 to 8 percent. A map of mean temperature along the line 
of sight can be computed from the ratio of the [0111] lines 5007 and 
4j63 after Usterbrock (Heference 10, p. 100) using the relation: 


TOO 


32900 


/i (5007 )\ 

l0ge \I(4363)/ 


This map is presented in Figure 23. 
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The electron density could be computed from the [SII] lines 6731 
and 6716, using a table calculated by Osterbrock (Reference 10, p. 112) 
for an assumed temperature of 10 4 K. This map is presented in Figure 24. 
The bright rings are artifacts. 

To illustrate the spatial distribution of material in different 
ionization states, the following three figures are presented. Figure 25 
is a mosaic of images in different bands, each color coded green, added 
to Ha color coded as red. The relative intensities have been adjusted 
to reveal the most information. In Figures 26 and 27, the ionization 
structure has been enhanced using the HSI transformation (Section VI) 
and assigning three bands to the color primaries blue, green, and red 
according to Table 2. The radial position of each ionization shell 
is apparent. 


Table 2. Color Assignments of Narrow Band Images for 
Figures 26 and 27 




Color Primaries 


Figure 


B1 ue 

Green 

Red 

26 

Hell, 4686 

Hel, 5876 
[QIII], 5007 

[01], 6300 

27 

[Oil], 3727 

[01], 6300 


These techniques should be applicable to any objects in an optically 
thin, non-LTE environment. 
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NGCt>720 mosaic of nine emission line images af'er radiometric 
correction and scaling for visual interpretation. Coherent , 
periodic, and random n'ise is clearly a problem in several 
of the images. 
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Figure 'b. NuCc>720 t wo-d imensional Fourier transforms displayed as 
log amplitude in a one-tu-one correspondence with the 
images of Figure 17 . The object data populates only the 
central region of each transform. Coherent noise is 
seen as clusters and periodic noise as streaks. 
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NGCP720 mosaic of modified transforms in the same format 
as Figure 16. The coherent noise spikes have been greatly 
reduced. Lach display has been renormalized to reveal 
noise residuals. 
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Figure 21. NGCt>720 mosaic 1’ polar coordinate transformations 

in the same order as Figure 20. Radial sections from 
the central star outward map as horizontal lines in this 
figure. A round nebula would appear as a vertical structure 
in this format. 
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Figure NuCb720 temperature map sealed foi viewing. The 3ky was 

tnresholded to zero because there the ratios of t tie luili] 
images at b007 A and A are undefined. 
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rigure 27. 


NGCb720 color-coded rendition of ionization structure. 
Color assignments are blue = [Oil], 3727; green = lOIII], 
5007; red = l 01], 6300. Colors have _een enhanced using' 
the HSI transformation. 
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SECTION VI 
COLOR ENHANCEMENTS 


In Sections IV and V of Reference 1, two powerful tools were 
illustrated for enhancing color differences between multispectral images 
of the same object. In this section, these two techniques (the HSI and 
Principal Component transformation) will be developed further and a third 
technique introduced, called the Log Intensity transformation. All of 
the images in this section were obtained from Dr. Allan Sandage of Hale 
Observatories. The images consisted of four spectral bands, 103a0 + WG2, 
103aD + GG14, 103aE + RG1, 103aE + Ha, of M82 obtained at the 200- inch 
telescope. The 103a0 image was synthesized from the average of four 
separate exposures, each registered and histogram matched. No correction 
was applied for the HD curve of the emulsion. Each image was linearly 
scaled so that the sky level was equal and the variance of the image 
within the region of the galaxy was also equal (this is equivalent to 
equating the volumes under the povvr spectra, thus forcing the information 
content of each band to be the same) . The latter constraint is certainly 
not a real one, but it provides a good basis for analysis because each 
image will be represented equally, none being dominated by brighter ones. 

In the remainder of this section, each of the three techniques will 
be applied to the same M82 images so that a comparative evaluation is 
possible . 


A. LOG INTENSITY STRETCH 

It is possible to vary the dynamic range of a set of multispectral 
images in any fashion whether it be linear or nonlinear and to preserve 
the color information (saturation and hue) by operating upon the intensity 
information only. In this application, it was convenient to reduce the 
dynamic range of the images such that they could be represented in 8 bits 
per pixel, a format compatible with the HSI and Principal Component 
transformations (and within the range of playback devices) to be addressed 
later. Because the majority of the galaxy resided at low intensities, a 
logarithmic compression was selected to reduce the upper dynamic range 
more than the lower dynamic range. Rather than take the logarithm of each 
pixel in turn (with disastrous consequences on the color) the logarithm 
was taken of the sum of each color band at each pixel, providing a scalar 
which could in turn be applied to that pixel location on each picture. 


! " 

*j = I j C i l0 6 c 2 2Z 



j=1,2 n (9) 


I. represents the intensity of band j of n bands at one pixel, and C 1 
and C 2 are constants which control the scaling and linearity respectively. 
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Figure 28 illustrates four oolor combinations of M82 obtained by select- 
ing the original four bands, three at a time. The color assignments are 
listed in Table 3, and details of the emulsion and filter combinations 
are shown in Table 4. 


Table 3» Standard Color Assignments for the M82 Mosaics 




Color Primaries 


Quadrant 

Blue 

Green 

Red 

Upper left 

103 a 0 

103aD 

103aE 

Upper right 

103a0 

103aD 

Ha 

Lower left 

103a0 

103aE 

Ha 

Lower right 

103aD 

103aE 

Ha 


Table 4. Etaulsion/Filter Combinations 




i 

i 


; \ 

'1 


Huulsion 

Filter 

103a0 

WG2 

103aD 

GG14 

103aE 

RG1 

103aE 

Ha 


Figure 28 consists of the original intensity images. One can 
distinguish much structure at the center of the galaxy, but not much 
is visible at lower intensities. Figure 29 illustrates the same mosaic 
as Figure 28, but after the Log Intensity transformation. Contrast 
is much lower along the bright bar but lower intensities are more readily 
observed. The color playback device can be operated under several cali- 
bration conditions. In the case of Figures 28 and 29, the device produced 
a print where equal steps in intensity in each band mapped into equal 
density steps in the photographic negative. This produces a pleasing 
dynamic range, segregates well between intensity levels in the print, but 
allows some cross talk between the colors so that hue becomes a function 
of intensity. Figure 30 is identical with Figure 29 except that the 
playback device was calibrated to produce a print where the reflectance 
of each step of the gray scale in each band (blue, green, or red) is 
proportional to the intensity of that band. This requirement provides 
an additional logarithmic compression effect, enhancing the lower intensity 
range. Colors should be truer in Figure 30 than in Figure 29. Both 
versions will also be presented later to aid in interpretation. Note 
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that the Log Intensity stretch has not produced a significant alteration 
of the colors in the scene; it acts as a means of redistributing intensity 
information only. Considerable filamentary detail can be seen in Figure 30, 
detail that remains invisible in Figure 28. 


B. HSI 

The HSI or Hue, Saturation, and Intensity transformation was 
introduced in References 1 and 11 where it was used to enhance color 
differences in KGC1097. Briefly, this transformation attempts to mimic 
the physiological parameterization of the color space from a polar 
coordinate transformation in which intensity is the sum of the three 
spectral bands (blue, green, and red), saturation is the colatitude 
from the gray axis, and hue is the longitude about the gray axis. In 
this section, the HSI transformation will be applied to M82 as a comparison 
with the Principal Component transformation in the next section, but 
also to illustrate some further manipulations of the HSI space. 

The images input to the HSI space are in this case the products of 
the Log Intensity transformation rather than the intensities themselves. 
Figure 31 illustrates a mosaic of the four spectral bands produced by the 
Log Intensity process, in the order upper left = 103a0, upper right = 103aD, 
lower left = 103aE, and lower right = Ha. Subsequent color products will 
be combinations taken three at a time as in Table 3 . The HSI processing 
consists of first generating the Hue, Saturation, and Intensity pictures 
from the blue, green, and red inputs. This is followed by performing a 
nonlinear contrast enhancement upon the Saturation picture so that the 
output histogram is equally populated everywhere. An optional step is 
introduced consisting of performing a low pass filter upon the histogram 
of the Hue picture and then nonlinearly contrast enhancing the Hue 
picture so that the Hue histogram would resemble the low pass filtered 
version. This last operation, although at first obscure, has a simple 
function. Concentrations of values in the Hue histogram are caused 
by regions of the pictures which have the same mean wavelength (the 
pattern recognition definition of an object "class"). It would be 
desirable to perform a higher th*»n normal contrast enhancement in these 
regions of the Hue histogram in order to improve the disoriminability 
of the eye to color differences within a single "class" of object which, 
because of the clustering of values, probably contains a great deal of 
information. This selective contrast enhancement is performed by first 
low pass filtering the Hue histogram, a process which smears or broadens 
the histogram clusters, and then using this histogram as a model for the 
desired output histogram. The contrast enhancement table generated to 
achieve the desired output histogram is highly nonlinear but is still 
monotonia so that crossovers of color are forbidden. Several advantages 
of this enhancement technique are that the filter preserves the means 
of the clusters; thus, the mean wavelength (mean hue) of the clusters 
is unaltered and the enhancement can be performed automatically. 

Figures 32 through 35 illustrate various stages in the above process 
for each of the four combinations of three colors selected as indicated 
in Table 5. 
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Table 5. Assignments of Bands to Primaries in Figures 32 through 35 


Figure 


Color Primaries 


Blue 

Green 

Red 

32 

103a0 

I03a0 

103aE 

33 

103a0 

103aD 

Ha 

34 

I03a0 

103aE 

Ha 

35 

103aD 

103aE 

Ha 


In each figure, the left column consists of three pictures in 
the following order: top is blue, middle is green, and bottom is red 
after the saturation and hue enhancements. These can be compared with 
the corresponding unenhanced versions in Figure 31. The right column 
in these same figures consists of enhanced HSI images in the following 
order: top is hue, middle is saturation, and bottom is intensity. 

The histograms at the bottom of the left column correspond one to one 
with the original unenhanced versions of the pictures above them. At 
the bottom of the right column of pictures are two columns of histograms. 
These are in order from top to bottom: hue, saturation, and intensity; 

the left column represents the original histograms before enhancement 
and the right column after enhancement. 

Spikes every 60° in the hue histogram are artifacts due to coinci- 
dences of planes of constant hue with the three-dimensional lattice 
of 8-bit color space. These spikes confuse the operation of the hue 
contrast enhancement and reduce its effectiveness. 

Color reconstructions based upon contrast enhancements of the 
saturation only (no modification of the hue) are presented in Figures 
36 and 37 in the order of Table 3. The modified hue images require 
careful interpretation and thus are not shown as color reconstructions. 
Figure 36 is calibrated with equal intensity steps corresponding with 
equal density steps and Figure 37 has equal intensity steps as equal 
reflectance steps as previously discussed. Considerable color structure 
is visible, betraying the spatial distribution of features with similar 
spectral properties. 

For some analyses, the intensity and saturation of a scene proves 
to be a nuisance. In order to see the mean wavelengths residing everywhere 
in the scene, the inverse HSI transformation was performed as usual except 
with both the intensity and saturation images discarded. This product 
is presented as Figure 38 in the usual Table 3 format. Some regions 
of the image which were too bright or dark to register color before 
are now visible. The sky is quite complex in hue because it has been 
raised above where the true sky level should be in order to observe 
faint structure, resulting in random noise of all spatial frequencies 


dominating the color. The hue of a pixel with either an intensity 
or a saturation of zero is undefined. 


C. PRINCIPAL COMPONENT 

The Principal Component transformation was introduced in Reference 
1 as a means of discriminating color structure without necessarily 
color coding it in a way traceable to the true color. A method has 
been developed (Reference 12) which allows a reconstruction of the 
approximate colors of a scene after enormous color exaggeration. The 
method does allow some cross talk between color channels and is thus 
dangerous to use if no care is taken in its analysis, but the technique 
provides the most powerful enhancement available. The reasons for 
this can be considered advantageous or disadvantageous depending upon 
one’s point of view. The reasons are as follows: 

(1) The algorithm performs equally well in discriminating color 
differences even if the images are not properly intensity 
balanced to begin with. 

(2) The technique is not limited to three colors, but can 
operate in a space of dimensionality N where N > 2. 

(3) The process is entirely automated, deriving the optimal 
transformation from diagonalizing the covariance matrix, 
a statistical quantity computed from the images. 

(4) The process formats images in order of descending information 
content, allowing manipulation of intermediate products. 

briefly, the Principal Component transformation consists of deter- 
mining that linear transformation (R) between the images which solves 
the set of matrix equations: 


R fc KR = K' (10) 

where K is the covariance matrix of the pictures, t means transpose, 
and K' is the matrix whose diagonal elements are the eigenvalue solutions 
to the characteristic matrix of K. The enhancement process consists of 
determining the transformation matrix (R) and then of transforming 
the different spectral bands with this matrix into uncorrelated images. 

A mosaic of the uncorrelated images is presented as Figure 39. 

The histograms below match one to one with the images above them. 
Considerable detail is evident in these images and it was at this stage 
that the application of this technique stopped in Reference 1. Table 
6 gives the eigenvalues for each of the images of M82 in Figure 39. 

The eigenvalues represent the variance or information content in each image. 
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Table 6. Information Content in Each Prinoipal Component 
Image of M82 for Figure 39 


Quadrant 

Eigenvalue 

Upper left 

4414 

Upper right 

679 

Lower left 

71 

Lower right 

33 


In order to understand the nature of the images of Figure 39 
better, the transformation matrix which generated them is presented. 


Upper left \ / 0.5118 0.5241 0.5070 0.4543X /103a0\ 

Upper right \ / -0.4578 -0.3456 0.0901 0.8141 \ / 103aD ) 

Lower left I s 1 -0.7136 0.4792 0.4469 -0.2474 / \103aE/(11) 

Lower right/ \ 0. 1377 -0.6132 0.7315 -0.2639/ \ Ha / 


For example, from the matrix, one can see that the upper right image 
is equivalent to Ha minus the average of 103a0 and 103aD, and that 
the upper left image is essentially the intensity transformation or 
sum of each image weighted by how well it correlates with the average. 

The color enhancement consists of contrast enhancing the principal 
component Images using a nonlinear enhancement which fills the output 
histogram so that it has a gaussian distribution, and then performing 
the inverse transformation using (R t ) back to a recognizable color 
space. In order to understand this process, it is important to visualize 
the color space. Because the original input images are highly correlated, 
the majority of the color •'pace (N-dimensional histogram) is populated 
along and around the graj ucis, leaving empty that portion of the apace 
far from the blue = green = . . . = red diagonal. Contrast enhancements 
upon the individual color bands assure that the entire diagonal is 
populated but it does little to populate off-axis regions of the color 
space. The Principal Component transformation performs a rotation 
of the color space such that what was once the gray axis (or the major 
information axis) now lies along one of the new axes, each new axis 
in turn being assigned an orthogonal direction through decreasingly 
smaller axes of the color space population until all degrees of freedom 
are exhausted. One can now see that the first principal component 
image is the intensity (upper left of Figure 39) because it coincides 
with the longest or intensity axis of the color space distribution. 

The narrowest axis through this distribution (lower right of Figure 39) 
is the least of all correlated directions and contains the least information 
and the most noise. Once the main axis of the color distribution has 
been aligned along a new coordinate direction, a contrast enhancement 
in the other coordinate direction (other pictures in Figure 39 than 
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th® upper left) fills the color spaoe. The inverse transformation 
using (R l ) re-rotates the oolor space to its original orientation such 
that blues return approximately to blues, eto. In the case of M82, 
all four bands were processed simultaneously in this fashion and the 
retransformed images were selected in groups of three for color addition 
in a manner similar to Table 3* Figure 40 illustrates the mosaic. 

Considerable color structure can be seen in the filaments. 

To reduce the noise level, each of the Principal Component images 
except the intensity image (upper left of Figure 39) was filtered using 
a two-dimensional median filter of a dimension which increased as the 
information content, given by the eigenvalue, decreased. Upon retransformation, 
the color-recombined products shown in Figure 41 appear softer and 
much less noisy. The equivalent operation performed upon the individual 
color bands would have soverly degraded the galaxy in order to reduce 
the noise an equivalent amount. 
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Hgure 28. M82 mosiac of four combinations of three colors ea^h. 

The images represent intensity information coded as equal 
density steps. All of the M82 mosaics appear in the order 
ot [able i in the text. The histograms beiow are of the 
sums of all three images. 







Figure 29. M62 mosiac after the Log Intensity transformation which preserves but 

does not enhance color differences. Images are coded so 'ha* log 
intensity is in equal aensity s'eps. See Tacit 3 f c * color designa- 
tions. The histograms are color coded for each spec’ral band assign- 
ment. The lower left histogram is * he sum of the ones above i*. 






Figure 30. M82 mosaic after Log Intensity transformation but coded 

with log intensity in equal reflectance steps. This implies 
a further logarithmic compression. See Table j for color 
designations. Filaments are readily visible in this version. 
Histograms are the same as for Figure 29- 
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Higure 3 1 * M62 mosaic of the four individual colors after the Leg 

Intensity transformation: (top left; I03a0 (blue); 

(top right) lOjali (green); (lower left) lQ3aE (broad 
band red); (lower right) lOjaE (Hi). The histograms 
are one to one with the images above. 






M82 mosaic of three spectral bands. The left column 
represents images after color enhancement from the hue 
saturation and intensity transformation. The right column 
represents top = hue, middle = saturation, bottom = intensity 
after contrast enhancement of the saturation and filtering 
of the hue histogram. Histograms beneath the right column 
of images represent hue, saturation, and intensity in the 
left column from top down and their enhancements in the 
right column. From top down, the left column pictures 
are I03a0, I03at>, and l03aE. 
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Figure 33- M8J mosaic in t. ho . e format as Figure 32. From the top 
down, the left column pictures are I03a0, itf3aD, and Hi. 
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M82 mosaic in the same format as Figure 32. From the top 
down, the left column pictures 3re I03a0, tQitb, and Hu. 
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Figure 3b. 


M82 mosaic after HS1 color enhancement with the hue, saturation, 
and intensity transformation and coded with log intensi f y in equal 
density steps. See Table 3 for color designations. The histograms 
below represent from left to right: (top row) hue, saturation, 

and intensity; (middle row) hue, contrast-enhanced saturation, and 
cont rast -enhanced intensity; (bottom row) blue, green, and red 
components of the picture. Compare this image with Figure 29. 
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M62 mosaic after HS1 color enhancement the same as figure 30 
except with log intensity coded into equal reflectance steps 
This version provides better filament detection and truer 
color consistency over the intensity range. See Table 3 
for color designations. Compare this image with Figure 30. 















Figure 39. 


L 


M82 mosaic of the principal component images produced from 
four spectral bands ILguation 11 ). The images are ranked 
in order of descending information content from upper left 
to lower right in rows as shown in Table t. Images differ 
because of color differences and, unline the original bands, 
are uncerrelated . 
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Figure 40. M82 mosaic of color-enhanced reconstructions from the principal 

component images of Figure 39- See Table 3 For color designa- 
t ons. Compare with Figures 36 and 37- 
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SECTION VII 


DECALIBRATION OP CCD ARRAYS 


An example of the procedure used to deoallbrate a CCD array la 
included in this section. The array in question was built by Texas 
Instruments as part of an investigation leading toward the development 
of the space telescope and Galileo sensors (Reference 13). The array 
suffers from many defects which will not be present in the final model 
and should not be considered as a perfected device. 

Figure 42 illustrates the raw image of a flat field produced by 
the CCD. The histograms illustrate the original intensity distribution 
at left and the contrast* enhanced distribution at right corresponding 
to the picture above. Several types of artifacts are present. The 
vertical white lines are due to damaged pixels which continually generate 
charge while the image is being read out vertically. The blaok column 
to the left is due to dead pixels which cannot transfer charge. The 
short vertical white lines at the top are due to leakage of charge back 
into the array from the transfer register at the top. The corners of 
the frame are saturated and are due to large blocks of pixels whioh 
continually generate charge (even in the dark current). Overall spatial 
structure is due to uneven thinning. The only important consideration 
is whether a pixel responds consistently to an exciting stimulus. If it 
does, it can be calibrated. Even the linearity of the response is of no 
importance; vidicons have been routinely calibrated at IPL (References 
14 and 15) when their response was nonlinear, spatially dependent, a 
function of the spatial structure of the scene and of previous scenes. 

In the case of the CCD, a first-order correction was applied 
because of the linearity of the device. No attempt was made to model 
the leakage to and from the transfer register, leaving this step to 
another algorithm not addressed here. The calibration consists of 
generating a large (greater than 15) number of flat field frames at 
known exposure values. Dark current frames are subtracted from each 
flat field. For each pixel, a straight line is fit by least squares 
to all of that pixel's values in all of the frames, relating exposure 
to response. The coefficients of this line are stored in a dataset 
for each pixel. Figures 43 and 44 show pictures of the slope and offset 
coefficients of a typical calibration file. To decalibrate a frame 
each pixel is reverse transformed using the appropriate coefficients 
and then scaled to the output. 

Figure 45 contains the decalibrated version of Figure 42. The 
histograms from left to right correspond to the original scene, the 
decalibrated picture scaled to the same mean value as the original frame, 
and the contrast-enhanced frame above. The histogram of the decalibrated 
frame is much narrower, approaching the ideal delta function of the true 
input stimulus. Only the random noise, integer truncation, and blemish 
structure in the image, which does not conform to the linear calibration 
model, give the histogram its final shape. The RMS error obtained from 
the output histogram agrees to within a few percent of the expected 
RMS error from theoretical considerations. The dark corners of the 


deoalibrated image and calibration files ara due to inadequate unsaturated 
data for fitting .'ie linear model. 

A display of the stability of each pixel is presented in Figure 46. 
This pioture oontains the varianoe of eaoh pixel as ooaputed from six flat 
field frames taken at identical exposure times. Large values indloate 
large unoertalnty (poor repeatability). Notloe that high varianoes do 
not necessarily correspond with blemishes in Figure 42. The large 
ourved oontours are due merely to the dynamic range available to the 
image. Certain pixels will appear less stable if the fractional value 
to be expressed in integer form lies dose to the deoision point midway 
between two integers. 

The use of leaat squares techniques in radiometric calibration 
can eliminate the noise from calibration files and dark currents alike, 
leaving the noise in the image as the only noise source. Another advantage 
of this approach la that it permits the use of higher order functions 
than linear to be used to model the light transfer function in the 
event that the device is nonlinear. This frees one from piecewise 
calibration files. 
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Figure i‘2. Scaled display of the response of a CCD array to a flat 
field stimulus. Large scale structure is due to uneven 
thinning. Vertical blemishes are due to dead pixels, 
continuously active pixels, and leakage from the transfer 
register into the picture. Hie jrny important consideration 
is repeatability. 
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Figure *4 j. Scaled display of the multiplicative term of the calibration 
file used to decal ibrate the CCD image of figure 4 j 
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Scaled display of the offset term of the calibration file 
used to decal ibrate the CCD image of Figure 42 
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Sealed decal ibrated CC1) flat field image. The histogram 
are from left to right (t) from the original flat field 
image (Figure 42), (2) from the deealibrated image, and 
13) from the cont rast-enhan ed deealibrated image above. 
Many pixels included in blemishes have also responded to 
the decalibration operation. 
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Figure 4o. Variance image displaying the variances of each pixel in 
a set of five identical exposures in an undecal ibrat ed 
field. Lower histograms represent the variance at left 
and the scaled variance of the above image at right. 
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SECTION VIII 


SUMMARY AND CONCLUSIONS 


From the preceding sections the following conclusions have been 
reached. 

(1) The Gull maximum entropy technique of restoring convolved 
imagery is not able to enhance resolution more than linear 
techniques. It may, however, prove valuable in analysis 
of objects whose true nature cannot be predicted because 
it generates images virtually free of artifacts. This is 
not the case for one-dimensional data. 

(2) Maps of the spatial distribution of polarization in galaxies 
can be generated from photographic data if the intensity 

of the images shows differences between intensities or sums 
of intensities that are greater than the RMS noise bars. 

(3) Temporal changes in objects observed over long-time baselines 
can be observed by differencing suitably normalized images. 

(H) Maps of the spatial distribution of ionization, temperature, 
and electron density can be generated for planetary nebulae. 

(5) Investigations of the three-dimensional nature of planetary 
nebulae are possible by observing the value and sign of 
the volume emission coefficient using Abel transformations 
or their equivalent. 

(6) The large dynamic range inherent in galaxies can be reduced 
into the narrower range available for display devices in 

such a way that colors in the original scene are not destroyed. 
A logarithmic intensity compression is advantageous for 
galaxies . 

(7) The Hue, Saturation, and Intensity transformation and the 
Principal Component transformation can both be used to 
enhance subtle color differences without scrambling the 
hues in the scene. The latter can optimize the visible 
information content of multispectral images by adjusting 
the transformation to the statistics of the scene. 

(8) CCD arrays can be modeled by least squares generated flat 
field calibration files which offer low noise storage, a 
selection of orders for the function to be used in fitting 
the light transfer curve, and freedom from piecewise-linear 
models. These files decalibrate images to theoretical 
device performance levels. 
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